---
title: "7-4 Days alive and free of hospital by day 28"
description: |
This outcome was defined as the number of days alive and free of hospital
from randomisation to 28 days, calculated as 28 minus the number of days
spent in hospital. All patients dying within 28 days will be assigned zero
free days.
author:
- name: James Totterdell
affiliation: University of Sydney
- name: Rob Mahar
affiliation: University of Melbourne
date: today
freeze: true
---
```{r}
#| label: pkgs
#| code-summary: Load packages
library (ASCOTr)
library (tidyverse)
library (lubridate)
library (kableExtra)
library (patchwork)
library (cmdstanr)
library (posterior)
library (bayesplot)
library (ggdist)
library (lme4)
library (broom)
library (broom.mixed)
library (bayestestR)
theme_set (theme_classic (base_size = 10 , base_family = "Palatino" ) +
theme (panel.grid = element_blank (),
strip.background = element_blank ()))
bayesplot_theme_set (theme_set (theme_classic (base_size = 10 , base_family = "Palatino" ) +
theme (panel.grid = element_blank (),
strip.background = element_blank ())))
color_scheme_set ("red" )
options (digits = 4 )
```
```{r}
#| label: analysis-sets
#| code-summary: Prepare datasets
all_data <- readRDS (file.path (ASCOT_DATA, "all_data.rds" )) |>
mutate (
out_dafh2 = if_else (D28_death == 0 & is.na (out_dafh), 28 - DD_total_days, out_dafh)
)
# FAS-ITT
fas_itt_dat <- ASCOTr::: make_fas_itt_set (all_data)
fas_itt_nona_dat <- fas_itt_dat |>
filter (! is.na (out_dafh))
# ACS-ITT
acs_itt_dat <- ASCOTr::: make_acs_itt_set (all_data)
acs_itt_nona_dat <- acs_itt_dat |>
filter (! is.na (out_dafh))
# AVS-ITT
avs_itt_dat <- ASCOTr::: make_avs_itt_set (all_data)
avs_itt_nona_dat <- avs_itt_dat |>
filter (! is.na (out_dafh))
```
```{r}
#| label: models
#| code-summary: Load models
ordmod0 <- compile_cmdstanr_mod (
file.path ("ordinal" , "logistic_cumulative" ), dir = "stan" )
ordmod <- compile_cmdstanr_mod (
file.path ("ordinal" , "logistic_cumulative_epoch_site" ), dir = "stan" )
ordmod_site <- compile_cmdstanr_mod (
file.path ("ordinal" , "logistic_cumulative_site" ), dir = "stan" )
logistic <- compile_cmdstanr_mod (
file.path ("binary" , "logistic_site_epoch" ), dir = "stan" )
```
```{r}
#| label: helper-functions
#| code-summary: Functions
make_summary_table_anticoagulation <- function (dat, format = "html" ) {
tdat <- dat %>%
group_by (CAssignment = factor (CAssignment,
levels = c ("C0" , "C1" , "C2" , "C3" , "C4" ),
labels = intervention_labels2 ()$ CAssignment)) %>%
summarise (
Patients = n (),
Known = sum (! is.na (out_dafh)),
Deaths = sprintf (
"%i (%.0f%%)" , sum (D28_death, na.rm = TRUE ), 100 * mean (D28_death, na.rm = TRUE )),
` DAFH, Median (Q1, Q3) ` = sprintf (
"%.0f (%.0f, %.0f)" ,
median (out_dafh, na.rm = T),
quantile (out_dafh, 0.25 , na.rm = TRUE ),
quantile (out_dafh, 0.75 , na.rm = TRUE ))
) %>%
bind_rows (
dat %>%
group_by (CAssignment = "Overall" ) %>%
summarise (
Patients = n (),
Known = sum (! is.na (out_dafh)),
Deaths = sprintf (
"%i (%.0f%%)" , sum (D28_death, na.rm = TRUE ), 100 * mean (D28_death, na.rm = TRUE )),
` DAFH, Median (Q1, Q3) ` = sprintf (
"%.0f (%.0f, %.0f)" ,
median (out_dafh, na.rm = T),
quantile (out_dafh, 0.25 , na.rm = TRUE ),
quantile (out_dafh, 0.75 , na.rm = TRUE ))
)
) %>%
rename (` Anticoagulation \n intervention ` = CAssignment)
kable (
tdat,
format = format,
align = "lrrrr" ,
booktabs = TRUE ,
linesep = ""
) %>%
kable_styling (
font_size = 9 ,
latex_options = "HOLD_position"
) %>%
row_spec (nrow (tdat), bold = T)
}
make_summary_table_antiviral <- function (dat, format = "html" ) {
tdat <- dat %>%
group_by (AAssignment = factor (AAssignment,
levels = c ("A0" , "A1" , "A2" ),
labels = intervention_labels2 ()$ AAssignment)) %>%
summarise (
Patients = n (),
Known = sum (! is.na (out_dafh)),
Deaths = sprintf (
"%i (%.0f%%)" , sum (D28_death, na.rm = TRUE ), 100 * mean (D28_death, na.rm = TRUE )),
` DAFH, Median (Q1, Q3) ` = sprintf (
"%.0f (%.0f, %.0f)" ,
median (out_dafh, na.rm = T),
quantile (out_dafh, 0.25 , na.rm = TRUE ),
quantile (out_dafh, 0.75 , na.rm = TRUE ))
) %>%
bind_rows (
dat %>%
group_by (AAssignment = "Overall" ) %>%
summarise (
Patients = n (),
Known = sum (! is.na (out_dafh)),
Deaths = sprintf (
"%i (%.0f%%)" , sum (D28_death, na.rm = TRUE ), 100 * mean (D28_death, na.rm = TRUE )),
` DAFH, Median (Q1, Q3) ` = sprintf (
"%.0f (%.0f, %.0f)" ,
median (out_dafh, na.rm = T),
quantile (out_dafh, 0.25 , na.rm = TRUE ),
quantile (out_dafh, 0.75 , na.rm = TRUE ))
)
) %>%
rename (` Antiviral \n intervention ` = AAssignment)
kable (
tdat,
format = format,
align = "lrrrr" ,
booktabs = TRUE ,
linesep = ""
) %>%
kable_styling (
font_size = 9 ,
latex_options = "HOLD_position"
) %>%
row_spec (nrow (tdat), bold = T)
}
```
# Outcome Derivation
The outcome is calculated for a patient as:
- missing, if day 28 mortality is unknown (`D28_PatientStatus` ) or if the patient was known to have been alive at day 28 but number of days in hospital (`D28_OutcomeTotalDaysHospitalised` ) is unknown
- 0 if they died by day 28 (`D28_PatientStatus` )
- 28 - min(28, `D28_OutcomeTotalDaysHospitalised` ) otherwise
As a cross check, @fig-dis-day-vs-days-hospitalised plots the number of days hospitalised (as reported in `D28_OutcomeTotalDaysHospitalised` ) against the day of discharge from the index admission.
There are some participants who were known to be alive at day 28, but had unknown total number of days hospitalised.
We may choose to assume that these patients were not re-admitted, and so their total number of days free of hospital is imputed as 28 - `DD_total_days` , that is, 28 minus the number of days spent in hospital during their index admission.
```{r}
#| label: fig-dis-day-vs-days-hospitalised
#| code-summary: Compare D28_OutcomeTotalDaysHospitalised to days hospitalised per dischage date.
#| fig-cap: |
#| Number of days hospitalised against day of discharge from index admission.
#| fig-height: 4
pdat <- fas_itt_dat %>%
filter (D28_death != 1 ) %>%
dplyr:: count (D28_day = pmin (28 , D28_OutcomeTotalDaysHospitalised), DIS_day = pmin (28 , DIS_day))
ggplot (pdat, aes (D28_day, DIS_day)) +
geom_point (aes (size = n), shape = 21 ) +
geom_abline () +
scale_x_continuous ("Number of days hospitalised (D28_OutcomeTotalDaysHospitalised)" ,
breaks = seq (0 , 30 , 2 )) +
scale_y_continuous ("Day of discharge from index admission" ,
breaks = seq (0 , 30 , 2 )) +
labs (size = "Count" )
```
# Descriptive {#descriptive}
The overall distribution of days alive and free of hospital (DAFH) are reported in @fig-dafh-dist for all participants in the AVS-ITT set.
```{r}
#| label: fig-dafh-dist
#| code-summary: Overall distribution of DAFH
#| fig-cap: |
#| Distribution of days alive and free of hospital, FAS-ITT.
#| fig-subcap: FAS-ITT
#| fig-height: 4
pdat <- fas_itt_nona_dat %>%
dplyr:: count (dafh = ordered (out_dafh, levels = 0 : 27 ), .drop = F) %>%
mutate (p = n / sum (n))
p <- ggplot (pdat, aes (dafh, p)) +
geom_col () +
labs (
x = "Days alive and free of hospital" ,
y = "Proportion"
)
path <- file.path ("outputs" , "figures" , "outcomes" , "secondary" )
ggsave (file.path (path, "outcome-7-4-dist-overall.pdf" ), p, height = 2.5 , width = 6 )
p
```
## Anticoagulation
```{r}
#| label: tbl-7-4-summary-anticoagulation
#| code-summary: Summary of DAFH outcome by arm
#| tbl-cap: Summary of days alive and free of hospital to day 28 by anticoagulation treatment group, FAS-ITT.
save_tex_table (
make_summary_table_anticoagulation (fas_itt_dat, "latex" ),
file.path ("outcomes" , "secondary" , "7-4-anticoagulation-summary" ))
make_summary_table_anticoagulation (fas_itt_dat)
```
```{r}
#| label: fig-dafh-dist-anticoagulation
#| code-summary: DAFH by intervention
#| fig-cap: |
#| Distribution of days alive and free of hospital by anti-coagulation intervention, FAS-ITT.
#| fig-height: 4
pdat <- fas_itt_nona_dat %>%
dplyr:: count (
CAssignment = factor (CAssignment, labels = str_replace (intervention_labels ()$ CAssignment, "<br>" , " \n " )),
dafh = ordered (out_dafh, levels = 0 : 27 ),
.drop = F) %>%
group_by (CAssignment) %>%
mutate (p = n / sum (n)) %>%
mutate (cp = cumsum (p)) %>%
ungroup ()
p <- ggplot (pdat, aes (CAssignment, p)) +
geom_col (aes (fill = fct_rev (dafh))) +
labs (x = "Anticoagulation intervention" , y = "Cumulative proportion" ) +
scale_fill_viridis_d ("Days alive and \n free of hospital" , option = "A" , direction = - 1 ) +
guides (fill = guide_legend (reverse = TRUE )) +
theme (legend.key.size = unit (0.75 , "lines" ))
pth <- file.path ("outputs" , "figures" , "outcomes" , "secondary" )
ggsave (file.path (pth, "outcome-7-4-descriptive-anticoagulation.pdf" ), p, height = 3 , width = 6 )
p
```
```{r}
#| label: fig-dafh-dist-anticoagulation2
#| code-summary: DAFH by intervention
#| fig-cap: |
#| Cumulative distribution of days alive and free of hospital by anti-coagulation intervention.
#| fig-height: 4
p <- ggplot (pdat, aes (dafh, cp)) +
geom_step (aes (colour = CAssignment, group = CAssignment)) +
labs (x = "Days alive and free of hospital" , y = "Cumulative proportion" ) +
scale_colour_viridis_d (
"Days alive and \n free of hospital" , option = "A" , begin = 0 , end = 0.8 ) +
guides (fill = guide_legend (reverse = TRUE ))
pth <- file.path ("outputs" , "figures" , "outcomes" , "secondary" )
ggsave (file.path (pth, "outcome-7-4-descriptive-anticoagulation2.pdf" ), p, height = 3 , width = 6 )
p
```
```{r}
#| label: fig-dafh-dist-anticoagulation-avs-itt
#| code-summary: DAFH by intervention
#| fig-cap: |
#| Distribution of days alive and free of hospital by anti-coagulation intervention, AVS-ITT.
#| fig-height: 4
pdat <- avs_itt_nona_dat %>%
dplyr:: count (
CAssignment = factor (CAssignment, labels = str_replace (intervention_labels ()$ CAssignment, "<br>" , " \n " )),
dafh = ordered (out_dafh, levels = 0 : 27 ),
.drop = F) %>%
group_by (CAssignment) %>%
mutate (p = n / sum (n)) %>%
mutate (cp = cumsum (p)) %>%
ungroup ()
p <- ggplot (pdat, aes (CAssignment, p)) +
geom_col (aes (fill = fct_rev (dafh))) +
labs (x = "Anticoagulation intervention" , y = "Cumulative proportion" ) +
scale_fill_viridis_d ("Days alive and \n free of hospital" , option = "A" , direction = - 1 ) +
guides (fill = guide_legend (reverse = TRUE )) +
theme (legend.key.size = unit (0.75 , "lines" ))
pth <- file.path ("outputs" , "figures" , "outcomes" , "secondary" )
ggsave (file.path (pth, "outcome-7-4-descriptive-anticoagulation-avs-itt.pdf" ), p, height = 3 , width = 6 )
p
```
## Antiviral
```{r}
#| label: tbl-7-4-summary-antiviral
#| code-summary: Summary of DAFH outcome by arm
#| tbl-cap: Summary of days alive and free of hospital to day 28 by antiviral treatment group, FAS-ITT.
save_tex_table (
make_summary_table_antiviral (fas_itt_dat, "latex" ),
file.path ("outcomes" , "secondary" , "7-4-antiviral-summary" ))
make_summary_table_antiviral (fas_itt_dat)
```
```{r}
#| label: fig-dafh-dist-antiviral
#| code-summary: DAFH by intervention
#| fig-cap: |
#| Distribution of days alive and free of hospital by antiviral intervention, FAS-ITT.
#| fig-height: 4
pdat <- fas_itt_nona_dat %>%
dplyr:: count (
AAssignment = factor (AAssignment, labels = str_replace (intervention_labels ()$ AAssignment, "<br>" , " \n " )),
dafh = ordered (out_dafh, levels = 0 : 27 ),
.drop = F) %>%
group_by (AAssignment) %>%
mutate (p = n / sum (n)) %>%
mutate (cp = cumsum (p)) %>%
ungroup ()
p <- ggplot (pdat, aes (AAssignment, p)) +
geom_col (aes (fill = fct_rev (dafh))) +
labs (x = "Antiviral intervention" , y = "Cumulative proportion" ) +
scale_fill_viridis_d ("Days alive and \n free of hospital" , option = "A" , direction = - 1 ) +
guides (fill = guide_legend (reverse = TRUE )) +
theme (legend.key.size = unit (0.75 , "lines" ))
pth <- file.path ("outputs" , "figures" , "outcomes" , "secondary" )
ggsave (file.path (pth, "outcome-7-4-descriptive-antiviral.pdf" ), p, height = 3 , width = 6 )
p
```
```{r}
#| label: fig-dafh-dist-antiviral2
#| code-summary: DAFH by intervention
#| fig-cap: |
#| Cumulative distribution of days alive and free of hospital by antiviral intervention.
#| fig-height: 4
p <- ggplot (pdat, aes (dafh, cp)) +
geom_step (aes (colour = AAssignment, group = AAssignment)) +
labs (x = "Days alive and free of hospital" , y = "Cumulative proportion" ) +
scale_colour_viridis_d (
"Days alive and \n free of hospital" , option = "A" , begin = 0 , end = 0.8 ) +
guides (fill = guide_legend (reverse = TRUE ))
pth <- file.path ("outputs" , "figures" , "outcomes" , "secondary" )
ggsave (file.path (pth, "outcome-7-4-descriptive-antiviral2.pdf" ), p, height = 3 , width = 6 )
p
```
```{r}
#| label: fig-dafh-dist-antiviral-avs-itt
#| code-summary: DAFH by intervention
#| fig-cap: |
#| Distribution of days alive and free of hospital by antiviral intervention, AVS-ITT.
#| fig-height: 4
pdat <- avs_itt_nona_dat %>%
dplyr:: count (
AAssignment = factor (AAssignment, labels = c ("Usual care" , "Nafamostat" )),
dafh = ordered (out_dafh, levels = 0 : 27 ),
.drop = F) %>%
group_by (AAssignment) %>%
mutate (p = n / sum (n)) %>%
mutate (cp = cumsum (p)) %>%
ungroup ()
p <- ggplot (pdat, aes (AAssignment, p)) +
geom_col (aes (fill = fct_rev (dafh))) +
labs (x = "Antiviral intervention" , y = "Cumulative proportion" ) +
scale_fill_viridis_d ("Days alive and \n free of hospital" , option = "A" , direction = - 1 ) +
guides (fill = guide_legend (reverse = TRUE )) +
theme (legend.key.size = unit (0.6 , "lines" ))
pth <- file.path ("outputs" , "figures" , "outcomes" , "secondary" )
fpth <- file.path (pth, "outcome-7-4-descriptive-antiviral-avs-tt.pdf" )
ggsave (fpth, p + coord_flip (), height = 2.5 , width = 6 )
system (sprintf ("pdftoppm %s %s -png" , fpth, gsub (".pdf" , "" , fpth)))
p
```
```{r}
#| label: fig-dafh-dist-antiviral-acs-itt
#| code-summary: DAFH by intervention
#| fig-cap: |
#| Distribution of days alive and free of hospital by antiviral intervention, ACS-ITT.
#| fig-height: 4
pdat <- acs_itt_nona_dat %>%
dplyr:: count (
AAssignment = factor (AAssignment, labels = str_replace (intervention_labels ()$ AAssignment, "<br>" , " \n " )),
dafh = ordered (out_dafh, levels = 0 : 27 ),
.drop = F) %>%
group_by (AAssignment) %>%
mutate (p = n / sum (n)) %>%
mutate (cp = cumsum (p)) %>%
ungroup ()
p <- ggplot (pdat, aes (AAssignment, p)) +
geom_col (aes (fill = fct_rev (dafh))) +
labs (x = "Antiviral intervention" , y = "Cumulative proportion" ) +
scale_fill_viridis_d ("Days alive and \n free of hospital" , option = "A" , direction = - 1 ) +
guides (fill = guide_legend (reverse = TRUE )) +
theme (legend.key.size = unit (0.75 , "lines" ))
pth <- file.path ("outputs" , "figures" , "outcomes" , "secondary" )
ggsave (file.path (pth, "outcome-7-4-descriptive-antiviral-acs-itt.pdf" ), p, height = 3 , width = 6 )
p
```
## Age
```{r}
#| label: fig-dafh-dist-age
#| code-summary: DAFH by age
#| fig-cap: |
#| Distribution of days alive and free of hospital by age group, FAS-ITT.
#| fig-height: 4
pdat <- fas_itt_nona_dat %>%
dplyr:: count (
agegrp = cut (AgeAtEntry, c (18 , seq (25 , 75 , 5 ), 100 ), include.lowest = T, right = F),
dafh = fct_rev (ordered (out_dafh, levels = 0 : 27 )),
.drop = F) %>%
group_by (agegrp) %>%
mutate (p = n / sum (n))
pdat2 <- pdat %>%
group_by (agegrp) %>%
summarise (n = sum (n))
p1 <- ggplot (pdat2, aes (agegrp, n)) +
geom_col (colour = "grey40" , fill = "grey40" ) +
geom_vline (xintercept = 60 , linetype = 2 ) +
labs (y = "Number of \n participants" ,
x = "Age at entry" ) +
geom_vline (xintercept = 8.5 , linetype = 2 ) +
theme (axis.text.x = element_text (hjust = 1 , angle = 45 ))
p2 <- ggplot (pdat, aes (agegrp, p)) +
geom_col (aes (fill = dafh)) +
labs (x = "Age" , y = "Cumulative \n proportion" ) +
scale_fill_viridis_d ("DAFH" , option = "A" , direction = - 1 ) +
guides (fill = guide_legend (reverse = TRUE )) +
geom_vline (xintercept = 8.5 , linetype = 2 ) +
theme (axis.text.x = element_text (hjust = 1 , angle = 45 )) +
theme (legend.key.size = unit (0.25 , "lines" ))
p <- p1 | p2
path <- file.path ("outputs" , "figures" , "outcomes" , "secondary" )
ggsave (file.path (path, "7-4-age.pdf" ), p, height = 2.5 , width = 6 )
p
```
## Sex
```{r}
#| label: fig-dafh-dist-sex
#| code-summary: DAFH by sex
#| fig-cap: |
#| Distribution of days alive and free of hospital by sex, FAS-ITT.
#| fig-height: 4
pdat <- fas_itt_nona_dat %>%
dplyr:: count (
Sex,
dafh = ordered (out_dafh, levels = 0 : 27 ),
.drop = F) %>%
group_by (Sex) %>%
mutate (p = n / sum (n)) %>%
mutate (cp = cumsum (p)) %>%
ungroup ()
pdat2 <- pdat %>%
group_by (Sex) %>%
summarise (n = sum (n))
p1 <- ggplot (pdat2, aes (Sex, n)) +
geom_col () +
labs (
y = "Number of \n participants" )
p2 <- ggplot (pdat, aes (Sex, p)) +
geom_col (aes (fill = fct_rev (dafh))) +
labs (x = "Sex" , y = "Cumulative \n proportion" ) +
scale_fill_viridis_d ("DAFH" , option = "A" , direction = - 1 ) +
guides (fill = guide_legend (reverse = TRUE )) +
theme (legend.key.size = unit (0.25 , "lines" ))
p <- p1 | p2
path <- file.path ("outputs" , "figures" , "outcomes" , "secondary" )
ggsave (file.path (path, "7-4-sex.pdf" ), p, height = 2.5 , width = 6 )
p
```
## Oxygen
```{r}
#| label: fig-dafh-dist-oxygen
#| code-summary: DAFH by oxygen
#| fig-cap: |
#| Distribution of days alive and free of hospital by oxygen requirement, FAS-ITT.
#| fig-height: 4
pdat <- fas_itt_nona_dat %>%
dplyr:: count (
supp_oxy2 = factor (supp_oxy2, labels = c ("Not required" , "Required" )),
dafh = ordered (out_dafh, levels = 0 : 27 ),
.drop = F) %>%
group_by (supp_oxy2) %>%
mutate (p = n / sum (n)) %>%
mutate (cp = cumsum (p)) %>%
ungroup ()
pdat2 <- pdat %>%
group_by (supp_oxy2) %>%
summarise (n = sum (n))
p1 <- ggplot (pdat2, aes (supp_oxy2, n)) +
geom_col () +
labs (
y = "Number of \n participants" ,
x = "Supplemental oxygen" )
p2 <- ggplot (pdat, aes (supp_oxy2, p)) +
geom_col (aes (fill = fct_rev (dafh))) +
labs (x = "Supplemental oxygen" , y = "Cumulative \n proportion" ) +
scale_fill_viridis_d ("DAFH" , option = "A" , direction = - 1 ) +
guides (fill = guide_legend (reverse = TRUE )) +
theme (legend.key.size = unit (0.25 , "lines" ))
p <- p1 | p2
path <- file.path ("outputs" , "figures" , "outcomes" , "secondary" )
ggsave (file.path (path, "7-4-oxygen.pdf" ), p, height = 2.5 , width = 6 )
p
```
## Country
```{r}
#| label: fig-dafh-dist-country
#| code-summary: DAFH by country
#| fig-cap: |
#| Distribution of days alive and free of hospital by country, FAS-ITT.
#| fig-height: 4
pdat <- fas_itt_nona_dat %>%
dplyr:: count (
Country = fct_infreq (Country),
dafh = ordered (out_dafh, levels = 0 : 27 ),
.drop = F) %>%
group_by (Country) %>%
mutate (p = n / sum (n)) %>%
mutate (cp = cumsum (p)) %>%
ungroup ()
pdat2 <- pdat %>%
group_by (Country) %>%
summarise (n = sum (n))
p1 <- ggplot (pdat2, aes (Country, n)) +
geom_col () +
labs (
y = "Number of \n participants" ,
x = "Country of enrolment" )
p2 <- ggplot (pdat, aes (Country, p)) +
geom_col (aes (fill = fct_rev (dafh))) +
labs (x = "Country" , y = "Cumulative \n proportion" ) +
scale_fill_viridis_d ("DAFH" , option = "A" , direction = - 1 ) +
guides (fill = guide_legend (reverse = TRUE )) +
theme (legend.key.size = unit (0.25 , "lines" ))
p <- p1 | p2
path <- file.path ("outputs" , "figures" , "outcomes" , "secondary" )
ggsave (file.path (path, "7-4-country.pdf" ), p, height = 2.5 , width = 6 )
p
```
## Site
```{r}
#| label: fig-dafh-dist-site
#| code-summary: DAFH by site
#| fig-cap: |
#| Distribution of days alive and free of hospital by study site, FAS-ITT.
#| fig-height: 4
pdat <- fas_itt_nona_dat %>%
dplyr:: count (
Country = factor (PT_CountryName, levels = c ("India" , "Australia" , "Nepal" , "New Zealand" ),
labels = c ("India" , "Australia" , "Nepal" , "New \n Zealand" )),
Site = fct_infreq (site),
dafh = ordered (out_dafh, levels = 0 : 27 )) %>%
complete (dafh = ordered (0 : 27 ), nesting (Country, Site), fill = list (n = 0 )) %>%
group_by (Country, Site) %>%
mutate (p = n / sum (n)) %>%
mutate (cp = cumsum (p)) %>%
ungroup () %>%
mutate (
Country = droplevels (Country),
Site = droplevels (Site)
)
pdat2 <- pdat %>%
group_by (Country, Site) %>%
summarise (n = sum (n)) %>%
ungroup ()
p1 <- ggplot (pdat2, aes (Site, n)) +
facet_grid ( ~ Country, scales = "free_x" , space = "free_x" ) +
geom_col () +
labs (
y = "Number of \n participants" ,
x = "" ) +
theme (axis.text.x = element_text (angle = 45 , hjust = 1 ),
panel.border = element_rect (fill = NA ))
p2 <- ggplot (pdat, aes (Site, p)) +
facet_grid ( ~ Country, scales = "free_x" , space = "free_x" ) +
geom_col (aes (fill = fct_rev (dafh))) +
labs (x = "Study Site" , y = "Cumulative \n proportion" ) +
scale_fill_viridis_d ("DAFH" , option = "A" , direction = - 1 ) +
guides (fill = guide_legend (reverse = TRUE , ncol = 1 )) +
theme (axis.text.x = element_text (hjust = 1 , angle = 45 )) +
theme (legend.key.size = unit (0.25 , "lines" ))
p <- p1 / p2 +
plot_layout (guides = 'collect' )
path <- file.path ("outputs" , "figures" , "outcomes" , "secondary" )
ggsave (file.path (path, "7-4-country-site.pdf" ), p, height = 4 , width = 6.25 )
p
```
## Calendar Time
```{r}
#| label: fig-dafh-dist-calendar
#| code-summary: DAFH by calendar date
#| fig-cap: |
#| Distribution of days alive and free of hospital by calendar time, FAS-ITT.
#| fig-height: 4
pdat <- fas_itt_nona_dat %>%
dplyr:: count (
yr = year (RandDate), mth = month (RandDate),
dafh = ordered (out_dafh, levels = 0 : 27 ),
.drop = F) %>%
group_by (yr, mth) %>%
mutate (p = n / sum (n)) %>%
mutate (cp = cumsum (p)) %>%
ungroup ()
p1 <- pdat %>%
group_by (yr, mth) %>%
summarise (n = sum (n)) %>%
ggplot (., aes (mth, n)) +
facet_grid ( ~ yr, drop = T, scales = "free_x" , space = "free" ) +
geom_col () +
labs (
y = "Number of \n participants" ,
x = "Calendar date (month of year)" ) +
scale_x_continuous (breaks = 1 : 12 )
p2 <- ggplot (pdat, aes (mth, p)) +
facet_grid ( ~ yr, drop = T, scales = "free_x" , space = "free" ) +
geom_col (aes (fill = fct_rev (dafh))) +
labs (x = "Calendar date (month of year)" , y = "Cumulative \n proportion" ) +
scale_fill_viridis_d ("DAFH" , option = "A" , direction = - 1 ) +
guides (fill = guide_legend (reverse = TRUE )) +
theme (legend.key.size = unit (0.25 , "lines" )) +
scale_x_continuous (breaks = 1 : 12 )
p <- p1 | p2
path <- file.path ("outputs" , "figures" , "outcomes" , "secondary" )
ggsave (file.path (path, "7-4-calendar-time.pdf" ), p, height = 2 , width = 6 )
p
```
## Sample Cumulative Logits
Some evidence of a shift from proportional odds at higher values of DAFH.
```{r}
#| code-summary: Plot sample cumulative logits
#| fig-cap: Inspect sample cumulative logits.
trt_counts <- fas_itt_nona_dat %>%
dplyr:: count (CAssignment, out_dafh) %>%
complete (CAssignment, out_dafh, fill = list (n = 0 )) %>%
group_by (CAssignment) %>%
mutate (p = n / sum (n))
trt_logit <- trt_counts %>%
group_by (CAssignment) %>%
mutate (clogit = logit (cumsum (p))) %>%
group_by (out_dafh) %>%
mutate (rel_clogit = clogit - mean (clogit)) %>%
filter (out_dafh != 27 )
ggplot (trt_logit, aes (out_dafh, rel_clogit)) +
facet_wrap ( ~ CAssignment) +
geom_point () +
labs (y = "Relative (to mean) sample cumulative logit" )
```
```{r}
#| code-summary: Plot sample cumulative logits
#| fig-cap: Inspect sample cumulative logits.
trt_counts <- fas_itt_nona_dat %>%
dplyr:: count (AAssignment, out_dafh) %>%
complete (AAssignment, out_dafh, fill = list (n = 0 )) %>%
group_by (AAssignment) %>%
mutate (p = n / sum (n))
trt_logit <- trt_counts %>%
group_by (AAssignment) %>%
mutate (clogit = logit (cumsum (p))) %>%
group_by (out_dafh) %>%
mutate (rel_clogit = clogit - mean (clogit)) %>%
filter (out_dafh != 27 )
ggplot (trt_logit, aes (out_dafh, rel_clogit)) +
facet_wrap ( ~ AAssignment) +
geom_point () +
labs (y = "Relative (to mean) sample cumulative logit" )
```
# Modelling
## FAS-ITT
The full platform model (excluding baseline CRP) is fit to the FAS-ITT dataset.
```{r}
#| label: fit-model-fas-itt
#| code-summary: Fit primary model
res <- fit_primary_model (
fas_itt_nona_dat,
ordmod,
outcome = "out_dafh" ,
intercept = FALSE
)
names (res$ drws$ AOR) <- "Nafamostat"
names (res$ drws$ COR) <- intervention_labels2 ()$ CAssignment[- (1 : 2 )]
names (res$ drws$ OR) <- c ("Ineligible aspirin" , "Age \u2265 60" , "Female" , "Oxygen requirement" , "Australia/New Zealand" , "Nepal" )
```
```{r}
#| label: odds-ratio-summary-table-fas-itt
#| code-summary: Odds ratio summary table
#| tbl-cap: Posterior summaries for model parameters (fixed-effects), FAS-ITT.
save_tex_table (
odds_ratio_summary_table_rev (c (res$ drws$ AOR, res$ drws$ COR, res$ drws$ OR), "latex" ),
"outcomes/secondary/7-4-primary-model-fas-itt-summary-table" )
odds_ratio_summary_table_rev (c (res$ drws$ AOR, res$ drws$ COR, res$ drws$ OR))
```
```{r}
#| label: fig-or-densities-fas-itt-model
#| fig-cap: Posterior densities for odds ratio contrasts.
p <- plot_or_densities (c (res$ drws$ AOR, res$ drws$ COR)) +
labs (x = "Odds ratio (log scale)" , y = "Comparison" )
pth <- file.path ("outputs" , "figures" , "outcomes" , "primary" , "7-4-primary-model-fas-itt-odds-ratio-densities.pdf" )
ggsave (pth, p, width = 6 , height = 2.5 )
p
```
```{r}
#| code-summary: Odds ratio summary for epoch and site
#| fig-cap: Summary of epoch and site posterior odds ratios.
p <- plot_epoch_site_terms (
res$ drws$ gamma_epoch,
res$ drws$ gamma_site,
factor (res$ dat$ region_by_site,
labels = c ("India" , "Australia \n New Zealand" , "Nepal" ))
)
pth <- file.path ("outputs" , "figures" , "outcomes" , "secondary" , "7-4-primary-model-epoch-site-terms-fas-itt.pdf" )
ggsave (pth, p, width = 6 , height = 4.5 )
p
```
:::{.callout-caution collapse="true"}
##### Diagnostics
```{r}
res$ fit$ summary (
c ("alpha" , "beta" , "gamma_epoch" , "gamma_site" , "tau_epoch" , "tau_site" )) %>%
print (n = Inf )
res$ fit$ diagnostic_summary ()
```
:::
:::{.callout-caution collapse="true"}
##### Trace plots
```{r}
mcmc_trace (res$ drws["beta" ])
mcmc_trace (res$ drws["alpha" ])
mcmc_trace (res$ drws["gamma_site" ])
mcmc_trace (res$ drws["gamma_epoch" ])
```
:::
### Posterior Predictive
```{r}
y_raw <- as.integer (levels (res$ dat$ y_raw))
y_ppc <- res$ drws$ y_ppc
y_ppc_raw <- rfun (\(x) y_raw[x])(y_ppc)
ppc_dat <- bind_cols (fas_itt_nona_dat, tibble (y_ppc = y_ppc_raw))
grp_ppc <- function (grp, d = 0 : 27 ) {
ppc_dat %>%
group_by (grp = {{grp}}) %>%
summarise (
y_lteq = map_dbl (d, ~ mean (out_dafh <= .x)),
ypp_lteq = map (d, ~ rvar_mean (y_ppc <= .x)),
y_eq = map_dbl (d, ~ mean (out_dafh == .x)),
ypp_eq = map (d, ~ rvar_mean (y_ppc == .x))
) %>%
unnest (c (ypp_lteq, ypp_eq)) %>%
mutate (days = d) %>%
ungroup () %>%
pivot_longer (
y_lteq: ypp_eq,
names_to = c ("response" , "event" ),
names_sep = "_" ,
values_to = "posterior" )
}
plot_grp_ppc <- function (dat, lab = "" , xlab = "Probability" ) {
ggplot (dat %>%
filter (response == "ypp" , event == "eq" ),
aes (y = days)) +
facet_wrap ( ~ grp, nrow = 1 , scales = "free_x" ) +
stat_slabinterval (aes (xdist = posterior), fatten_point = 1 ) +
geom_point (data = dat %>% filter (response == "y" , event == "eq" ),
aes (y = days, x = mean (posterior)),
colour = "red" ,
shape = 23 ) +
scale_x_continuous (
xlab, breaks = c (0 , 0.3 , 0.6 ),
sec.axis = sec_axis (
~ . , name = lab, breaks = NULL , labels = NULL )) +
scale_y_continuous ("Days" ) +
theme (strip.text = element_text (size = rel (0.7 )),
axis.title.x = element_text (size = rel (0.7 )),
axis.text.x = element_text (size = rel (0.65 )),
axis.title.y = element_text (size = rel (0.75 )),
axis.title.x.bottom = element_blank ())
}
pp_A <- grp_ppc (AAssignment)
pp_C <- grp_ppc (CAssignment)
pp_ctry <- grp_ppc (Country)
pp_epoch <- grp_ppc (epoch)
pp_site <- grp_ppc (site) %>%
left_join (ppc_dat %>% dplyr:: count (site, Country),
by = c ("grp" = "site" ))
p0 <- plot_grp_ppc (pp_A, "Antiviral" , "" )
p1 <- plot_grp_ppc (pp_C, "Anticoagulation" , "" )
p2 <- plot_grp_ppc (pp_ctry, "Country" , "" )
p3 <- plot_grp_ppc (pp_epoch, "Epoch" , "" )
p4 <- plot_grp_ppc (
pp_site %>% filter (Country == "IN" ), "Sites India" , "" )
p5 <- plot_grp_ppc (
pp_site %>% filter (Country == "AU" ), "Sites Australia" , "" )
p6 <- plot_grp_ppc (
pp_site %>% filter (Country == "NP" ), "Sites Nepal" , "" )
p7 <- plot_grp_ppc (
pp_site %>% filter (Country == "NZ" ), "Sites New Zealand" , "" )
g1 <- (p0 | p1 | p2) / p3
g2 <- p4 / p5 / (p6 | p7)
pth1 <- file.path (
"outputs" , "figures" , "outcomes" , "secondary" ,
"7-4-primary-model-fas-itt-ppc1.pdf" )
pth2 <- file.path (
"outputs" , "figures" , "outcomes" , "secondary" ,
"7-4-primary-model-fas-itt-ppc2.pdf" )
ggsave (pth1, g1, width = 6 , height = 5 , device = cairo_pdf)
ggsave (pth2, g2, width = 6 , height = 6 , device = cairo_pdf)
g1
g2
```
## ACS-ITT
```{r}
#| label: fit-model-acs-itt
#| code-summary: Fit primary model - ACS-ITT
res <- fit_primary_model (
acs_itt_nona_dat,
ordmod,
outcome = "out_dafh" ,
intercept = FALSE
)
names (res$ drws$ AOR) <- "Nafamostat"
names (res$ drws$ COR) <- intervention_labels2 ()$ CAssignment[- (1 : 2 )]
names (res$ drws$ OR) <- c ("Ineligible aspirin" , "Age \u2265 60" , "Female" , "Oxygen requirement" , "Australia/New Zealand" , "Nepal" )
```
```{r}
#| label: odds-ratio-summary-table-acs-itt
#| code-summary: Odds ratio summary table
#| tbl-cap: Posterior summaries for model parameters (fixed-effects), ACS-ITT.
save_tex_table (
odds_ratio_summary_table_rev (c (res$ drws$ AOR, res$ drws$ COR, res$ drws$ OR), "latex" ),
"outcomes/secondary/7-4-primary-model-acs-itt-summary-table" )
odds_ratio_summary_table_rev (c (res$ drws$ AOR, res$ drws$ COR, res$ drws$ OR))
```
```{r}
#| label: fig-or-densities-acs-itt-model
#| fig-cap: Posterior densities for odds ratio contrasts, ACS-ITT.
p <- plot_or_densities (c (res$ drws$ AOR, res$ drws$ COR)) +
labs (x = "Odds ratio (log scale)" , y = "Comparison" )
pth <- file.path ("outputs" , "figures" , "outcomes" , "primary" , "7-4-primary-model-acs-itt-odds-ratio-densities.pdf" )
ggsave (pth, p, width = 6 , height = 2.5 )
p
```
## AVS-ITT
### Pre-Specified Model
- excluding epoch
```{r}
#| label: fit-model-avs-itt-pre-spec
#| code-summary: Fit primary model - AVS-ITT
res <- fit_primary_model (
avs_itt_nona_dat |> mutate (ctry = droplevels (ctry)),
ordmod_site,
outcome = "out_dafh" ,
vars = c ("agegte60" , "sexF" , "supp_oxy2" , "crp_tertile" , "ctry" ),
beta_sd_var = c (2.5 , 2.5 , 2.5 , 2.5 , 2.5 , 2.5 , 1 ),
intercept = FALSE
)
names (res$ drws$ AOR) <- "Nafamostat"
names (res$ drws$ COR) <- intervention_labels2 ()$ CAssignment[- (1 : 2 )]
names (res$ drws$ OR) <- c ("Age \u2265 60" , "Female" , "Required oxygen" , "CRP (2nd tertile)" , "CRP (3rd tertile)" , "CRP (unknown)" , "Nepal" )
```
```{r}
#| label: odds-ratio-summary-table-avs-itt-pre-spec
#| code-summary: Odds ratio summary table
#| tbl-cap: Posterior summaries for model parameters (fixed-effects), AVS-ITT.
save_tex_table (
odds_ratio_summary_table_rev (c (res$ drws$ AOR, res$ drws$ COR, res$ drws$ OR), "latex" ),
"outcomes/secondary/7-4-primary-model-avs-itt-summary-table-pre-spec" )
odds_ratio_summary_table_rev (c (res$ drws$ AOR, res$ drws$ COR, res$ drws$ OR))
```
```{r}
#| label: fig-or-densities-avs-itt-model-pre-spec
#| fig-cap: Posterior densities for odds ratio contrasts, AVS-ITT.
p <- plot_or_densities (c (res$ drws$ AOR)) +
labs (x = "Odds ratio (log scale)" , y = "Comparison" )
pth <- file.path ("outputs" , "figures" , "outcomes" , "primary" , "7-4-primary-model-avs-itt-odds-ratio-densities-pre-spec.pdf" )
ggsave (pth, p, width = 6 , height = 2.5 )
p
```
```{r}
#| code-summary: Odds ratio summary for epoch and site
#| fig-cap: Summary of epoch and site posterior odds ratios.
p <- plot_site_terms (
res$ drws$ gamma_site,
factor (res$ dat$ region_by_site,
labels = c ("Australia \n New Zealand" , "Nepal" ))
)
pth <- file.path ("outputs" , "figures" , "outcomes" , "secondary" , "7-4-primary-model-site-terms-avs-itt-pre-spec.pdf" )
ggsave (pth, p, width = 6 , height = 4.5 )
p
```
### Reduced Model
- excluding epoch, site, and country due to small sample size
```{r}
#| label: fit-model-avs-itt
#| code-summary: Fit primary model - AVS-ITT
res <- fit_primary_model (
avs_itt_nona_dat,
ordmod0,
outcome = "out_dafh" ,
vars = c ("agegte60" , "sexF" , "supp_oxy2" , "crp_tertile" ),
beta_sd_var = c (2.5 , 2.5 , 2.5 , 2.5 , 2.5 , 2.5 ),
intercept = FALSE
)
names (res$ drws$ AOR) <- "Nafamostat"
names (res$ drws$ COR) <- intervention_labels2 ()$ CAssignment[- (1 : 2 )]
names (res$ drws$ OR) <- c ("Age \u2265 60" , "Female" , "Required oxygen" , "CRP (2nd tertile)" , "CRP (3rd tertile)" , "CRP (unknown)" )
```
```{r}
#| label: odds-ratio-summary-table-avs-itt
#| code-summary: Odds ratio summary table
#| tbl-cap: Posterior summaries for model parameters (fixed-effects), AVS-ITT.
save_tex_table (
odds_ratio_summary_table_rev (c (res$ drws$ AOR, res$ drws$ COR, res$ drws$ OR), "latex" ),
"outcomes/secondary/7-4-primary-model-avs-itt-summary-table" )
odds_ratio_summary_table_rev (c (res$ drws$ AOR, res$ drws$ COR, res$ drws$ OR))
```
```{r}
#| label: fig-or-densities-avs-itt-model
#| fig-cap: Posterior densities for odds ratio contrasts, AVS-ITT.
p <- plot_or_densities (c (res$ drws$ AOR)) +
labs (x = "Odds ratio (log scale)" , y = "Comparison" )
pth <- file.path ("outputs" , "figures" , "outcomes" , "primary" , "7-4-primary-model-avs-itt-odds-ratio-densities.pdf" )
ggsave (pth, p, width = 6 , height = 2.5 )
p
```
<!-- ### Re-analyse using imputed when outcome was missing -->
<!-- Impute missing values assuming no readmissions occurred. -->
<!-- ```{r} -->
<!-- #| label: fit-model-avs-itt-impute -->
<!-- #| code-summary: Fit primary model - AVS-ITT -->
<!-- res <- fit_primary_model( -->
<!-- avs_itt_dat |> filter(!is.na(out_dafh2)), -->
<!-- ordmod0, -->
<!-- outcome = "out_dafh2", -->
<!-- vars = c("agegte60", "sexF", "supp_oxy2", "crp_tertile"), -->
<!-- beta_sd_var = c(2.5, 2.5, 2.5, 2.5, 2.5, 2.5), -->
<!-- intercept = FALSE -->
<!-- ) -->
<!-- names(res$drws$AOR) <- "Nafamostat" -->
<!-- names(res$drws$COR) <- intervention_labels2()$CAssignment[-(1:2)] -->
<!-- names(res$drws$OR) <- c("Age \u2265 60", "Female", "Required oxygen", "CRP (2nd tertile)", "CRP (3rd tertile)", "CRP (unknown)") -->
<!-- ``` -->
<!-- ```{r} -->
<!-- #| label: odds-ratio-summary-table-avs-itt-impute -->
<!-- #| code-summary: Odds ratio summary table -->
<!-- #| tbl-cap: Posterior summaries for model parameters (fixed-effects), AVS-ITT. -->
<!-- odds_ratio_summary_table_rev(c(res$drws$AOR, res$drws$COR, res$drws$OR)) -->
<!-- ``` -->
<!-- ```{r} -->
<!-- #| label: fig-or-densities-avs-itt-model-impute -->
<!-- #| fig-cap: Posterior densities for odds ratio contrasts, AVS-ITT. -->
<!-- p <- plot_or_densities(c(res$drws$AOR)) + -->
<!-- labs(x = "Odds ratio (log scale)", y = "Comparison") -->
<!-- p -->
<!-- ``` -->
### Treatment-only
```{r}
#| label: fit-model-avs-itt-trt-only
#| code-summary: Fit primary model - AVS-ITT
res <- fit_primary_model (
avs_itt_nona_dat,
ordmod0,
outcome = "out_dafh" ,
vars = NULL ,
beta_sd_var = NULL ,
intercept = FALSE ,
includeC = FALSE
)
names (res$ drws$ AOR) <- "Nafamostat"
```
```{r}
#| label: odds-ratio-summary-table-avs-itt-trt-only
#| code-summary: Odds ratio summary table
#| tbl-cap: Posterior summaries for model parameters (fixed-effects), AVS-ITT.
odds_ratio_summary_table_rev (c (res$ drws$ AOR))
```
```{r}
#| label: fig-or-densities-avs-itt-model-trt-only
#| fig-cap: Posterior densities for odds ratio contrasts, AVS-ITT.
p <- plot_or_densities (c (res$ drws$ AOR)) +
labs (x = "Odds ratio (log scale)" , y = "Comparison" )
p
```